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SlWtARY 


The direct simulation methods developed by Orszag and Patterson (1972) for 
isotropic turbulence have been extended to homogeneous turbulence in an 
incompressible fluid subjected to uniform deformation or rotation. The results of 
simulations for irrotational strain (plane and axisymmetric ) » shear, rotation, and 
relaxation toward Isotropy following axisymmetric strain are compared with linear 
theory and experimental data. Emphasis is placed on the shear flow because of its 
importance and because of the availability of accurate and detailed experimental 
data. The computed results arc used to assess the accuracy of two popular models 
used in the closure of the Reynolds-stress equations. Data from a variety of the 
computed fields and the details of the numerical methods used in the simulation are 
also presented* 


INTRODUCTION 


The turbulence problem has renvained tlie greatest challenge to fluid dynamicists 
since its definition by Reynolds 100 years ago. The nonlinearity of the 
Navler-Stokes equations prevents statistical analysis of the dynamics of energetic 
turbulent fields, and as a result no theory has yet been devised that is devoid of 
assumptions beyond that of the validity of the Navicr-Stokes equations themselves. 
Although statistical, censorial, and dimensional analysis combined with the 
continuity condition provide constraints on the statistics of a turbulent field, our 
knowledge ot the statistics themselves has been gained primarily through experiment. 
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The method of direct numerical simulation of homogeneous turbulence in general 
use today Is that of Orszag and Patterson (1972). The turbulence field Is 
represented by the coefficients of a truncated three-dimensional Fourier scries, or 



equlvAlencIy by spAtially periodic discrete values on a uniform mesh in physical 
space. Trans f onset ion between the two spaces (physical and wave) Is accomplished by 
fast Fourier transforms. Differentiation Is performed in wave space as 
multiplication by the wave number » and velocity products are carried out in physical 
space to avoid the expensive convolution sums required for that operation in wave 
space. The field can be tlme**advanced In either space by any algorithm appropriate 
for systems of ordinary differential equations. 

There are two basic requirements that a direct simulation must meet If It is to 
represent turbulence. First, It must represent a solution of the Navler-Stokes 
equations. This means that all scales of motion must be adequately resolved, In the 
deterministic sense, by the computational mesh. In particular, the Reynolds number 
should be small enough to allow the mesh to accurately capture the viscous 
dissipation scales. The explicit use of Fourier series then Implies a spatially 
periodic solution of the Navler^^Stokes equation. The second requirement Is that this 
periodic solution be sufficiently complicated, that Is, It must provide adequate 
statistical resolution (large sample) of the set of all possible fluid motions 
allowed by the Navler*>Stokes equations. The computational sample of a scale of 
motion is Inversely proportional to the volume of the scale. Thus scales of motion 
comparable to the computational period have a very small sample (^1), and scales 
comparable to the mesh spacing have a large sample (^10^ for a 128^ mesh). The two 
requirements for a turbulence simulation conflict j the sample improves as the energy 
moves to smaller scales but the viscous resolution is degraded. 

The range of scales (the ratio of the wavelength of the fundamental to that of 
the highest wave number carried) in a given direction is M/2, where M is the number 
of nodes in that direction. The largest computers available today allow 128 mesh 
cells in each direction, which limits the useful range of scales to about 10 or so. 
That is, the bulk of the total turbulent energy must be contained within one decade 
of scale. This is a severe constraint and it limits completely resolved direct 
simulation to very low Reynolds numbers. In practice we usually accept some error in 
order to obtain a higher Reynolds number. In some cases we sacrifice sample and 
shift Che computed scale range toward the small scales , and in others we sacrifice 
resolution of the dissipation process and shift the computed range toward the large 
scales. For the current study we are most Interested in the anisotropy of the 
Reynolds-stress tensor, and less interested in the details (such as Intermittency ) of 
the small scales, and we generally accept incomplete resolution of the dissipation 
scales. It appears, however, that the error in total dissipation is much less 
affected than is its spectral distribution because the total is determined primarily 
by energy transfer down*-scale within the energetic scales; this in turn depends on 
how well the energetic scales themselves are resolved and not on how well the actual 
dissipative scales are resolved. 

Computation and physical experiment complement each other rather well for 
homogeneous turbulence. The experimental difficulties are associated with setting up 
the mean motion, achieving a homogeneous turbulence field, and in the actual 
measurement of the various spectra, correlations, etc. There is no doubt that the 
physical experiment produces results at all scales In accordance with the 
Navier-Stokes equations, although frequency response (in both space and time) of the 
Instrumentation can limit their measurement* The computation on the other hand has 
no difficulty setting up the mean flow or making measurements. The difficulty is 
that because of limited resolution (statistical accuracy at the large scales, 
numerical accuracy at the small scales), it is not always simple to relate computed 
quantities to those of the experiment at much higher Reynolds number. The goal of 
this study Is to use experimental data to validate simulated turbulence fields, which 
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In turn provide meaeureraente that are not available from experiment. The q^ntltlea 
measured here are primarily those volume averages appearing In the Reynolds-stress 

equations. 


In the sections that follow we consider four cases of homogeneous turbulence, 
each being the evolution from an Initially Isotropic state caused by a mean velocity 
gradient uniform In space and time. The deformations considered are (1) Irrotatlonal 
plane strain, (2) irrotatlonal axisymmetrlc strain, (3) uniform shear, and (4) 
uniform rotation. In the case of axisymmetrlc strain Ve allow the turbulence to 
relax back toward Isotropy when the straining ceases. The only theory available to 
describe the evolution of these flows is a linear theory valid in the limit in which 
the time scale of all turbulence scales la much longer than the time scale given by 
the mean velocity gradient. This theory, due originally to Taylor (1^35), has been 
developed for each case by various authors and is used, togethar with available 
experimental data, to establish (we hope) some credibility for the simulations. 
Then, as an example of how the detailed measurements of the sicvulated fields might be 
used to aid the turbulence-model maker, we have compared measured values from the 
fields with modeled values following recommendations of Rotta (1951) and Launder, 
Reece, and Rodl (1975)* 


Finally, more details are given in the appendixes. Appendix A presents in 
tabular form various raw data measured in each field, including the Reynolds-stress 
budget; Appendix B provides a more detailed exposition of the equations and numerical 
methods used in the simulations. 


PLANE STRAIN 


Nearly isotropic turbulence, subjected to a uniform strain in two directions, 
has been investigated experimentally by Townsend (1951) and by Tucker and Reynolds 
(1968) bv passing grid-generated turbulence through a channel of changing cross 
section.' Townsend imposed a total strain of 4 in directions transverse to the 
stream; Tucker and Reynolds imposed total strains of 6 in two different 
configurations, the first, like that Townsend, was strained In the transvers^ 

directions and the second had one strain Imposed in the flow direction. We follow 
the nomenclature of Townsend and take the flow direction to be x, imposing positive 
strain in the z direction and negative strain in the y direction. The results of 

Tucker and Reynolds Indicate no significant difference between the two types of 

strain. The longitudinal strain does have the advantage, like the axisymmetrlc 
strain discussed in the next section, that the strain is given by direct measurement 

of the variation of mean flow velocity down the channel. 

Results of four computed cases, all evolving from the same isotropic initial 
state are presented. The initial state itself is the result of an Isotropic 
simulation so that the spectral transfer is established. The microscale Reynolds 
number q\ii/v of the Initial field is *35, and dimensionless stra^rates 
12 and 4 vwhere a - “ -dv/dy, q2 - uiui, e - 2vui,Jui,j) 

are Imposed In cases B2D1 , B2D2. B2D3, and B2D4, respectively. The total strain 
achieved is 4; higher strains result in Inadequate resolution because tne 
computational mesh moves with the mean flow and becomes excessively strained. At the 
initial state, the computational cell has Av ■ 2Az, and at a strain of , z “ y. 
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Three Independent dimensionless parameters can be formed from the time t, strain 
rate a viscosity v, and turbulent velocity and length scales q and L respectively. 
These are taken to be at, aL/.q, qL/v, the dimensionless time, ratio of 
turbulence and strain time scales, and a turbulence macroscale Reynolds number. No 
information on macroscales is directly available from the experiments, so we have 
instead estimated L by assuming e -q^/L in the Isotropic initial state. The initial 
conditions are then scaled by the dimensionless groups aq^/e and q /ev. Estimates 
of the quantities for the experiments are as follows; 



aq^/ e 

e‘*/ve 

Townsend .5 in. 

mesh .45 

95 

1 in. 

mesh .45 

230 

Tucker-Reynolds 

.3 

425 

Present results 

.5, 1,2, 4 

129 


Here the isotropic relation e *“ lOvq^/Ajj has been used. For the Townsend 
experiment, the microscales are given and these are used to estimate for the 
Tucker and Reynolds experiment we have estimated e by fitting a power raw, with 
exponent of -1.2, to the plotted energy history; these estimates are very rough 
Indeed. 

The linear theory of "sudden distortion" was originated by Taylor (1935) who 
considered the invlscid irrotational strain of a typical Fourier component of the 
velocity field. The effect of viscosity can be Included by use of an Integrating 
factor (Pearson (1959)). Batchelor and Proudman (1954) extended the Invlscid theory 
to determine the effect of rapid irrotational strains on the energy tensor of an 
initially isotropic turbulence field, and the result is Independent of the form of 
the energy spectrum E(k). Although viscosity can in principle be : -ndled exactly , 
the necessary integrals over Fourier space contain exponential factors and probably 
cannot be evaluated in closed form. The result would depend on the initial spectrum. 

A comparison of the linear theory, experimental, and computational results is 
presented in figures 1 through 3. The evolution of the normalized energy tensor is 
shown in figure 1, and the structure parameters Introduced by Townsend are shown in 
figures 2 and 3. The computational results approach the linear theory consistently 
with increasing strain rate but, although they follow the trends of the experiment, 
quantitative agreement (particularly for Kj and Kf ) is lacking. 

The differences between the two experiments are not consistent with our 
estimates of their dimensionless time scales and Reynolds numbers. At a given level 
of strain, the structure parameter Kj should increase with strain rate but decrease 
wl^th Increasing Reynolds number. Our estimates of the dimensionless strain rates and 
Reynolds numbers of the experiments would cause us to expect higher Ki values from 
the the Townsend experiment than from the Tucker and Reynolds experiment, but that is 
not the case. 
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AXISYWETRIC STRAIN 


Turbulence in an axisymmetric strain field Is the simplest anisotropic flow. 
Because it occurs often in engineering appllcatio,)s, it has been studied ir a number 
of experiments. We will compare the computational results with the experimental data 
of Uberoi (1956) and Hussain and Ramjee (1976), who achieved total strains up to 16 
at a wide range of Reynolds numbers. Unfortunately, Information concerning the 
dissipation rate and length scales was n6t presented, and we are unable to determine 
with any confidence the turbulence Reynolds numbers or time scales of the turbulence 
prior to straining. 

Simulations were made with four constant strain rates applied at each of two 
viscosity values. The dimensionless strain rates aq^/e, where 
a - du/dx - -2dv/dy = -2dw/dz, were 0.3, 0.6, 1.2, and 2.4, and the initial 
Reynolds numbers q‘*/vc were 15 and 56. The evolution of the longitudinal and 
transverse component energies, normalized by their initial values for five of the 
eight simulated cases, the inviscld linear theory, and the experlmantal results of 
Uberoi and Hussaln-Ramjee , are shown in figure 4. The computational results approach 
rapid-distortion theory consistently with increasing strain rate and Reynolds number, 
and the differences between computation and experiment seem to indicate a higher 
experlmei^al strain rate. Beyond a strain of 4, the experimental component 
energy u2 grows because of the redistribution of energy by the pressure-strain 
correlation, and the computational data (appendix A) indicate that this is a result 
primarily of the growth of the "slow-pressure-strain" term with increasing strain. 


HOMOGENEOUS SHEAR 


Homogeneous shear turbulence is the closest structurally to flows of engineering 
interest and has been extensively investigated by experimenters. In a series of 

Hopkins University, Rose (1966), Champagne, Harris and Corrsln 
(1970), Harris, Graham and Corrsin (1977), and Tavoularis and Corrsln (1981) 

(hereafter TC) have provided extensive documentation of the development of turbulence 
subjected to uniform shear. The results of Mulhearn and Luxton (1975) are in 

substantial agreement with the earlier work of the Johns Hopkins group. On the other 

hand, the linear theory of horogeneous shear turbulence has been worked out by 

Deissler (1961, 1970, 1972), and Fox (1964); their results are consistent with the 
Initial stage of development of the experimentally observed turbulence. 

The linear theory can be worked out in a somewhat simpler way than that of 
Deissler, who worked directly with the equations for the spectrum tensor after 
dropping out the triple correlations. Written in a coordinate system moving with the 
mean shear, the Unearlzod Navier-Stokes equations are 

SV -f pjj * 

^ Py “ V'i|)v 

Wj. + pjp * vSPw 

Ux “ Vy -f W 2 * 0 
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Upon taking the Fourier transform of these equations and solving the remaining 
ordinary differential equations for the time variation, we obtain 


Fu ■ “Akj + sjks^ tan"^n - ki^ 1 

Fv - Bki 

Fv> - Ak^ - Bk^k3 |tan"^n + r—T — 2I 


where F is th« integrating factor 


log . i 

k2 - k^st 

n * 

^72 + ks2 


and A(k), B(k) are determined from the initial conditions. The evolution of the 
spectrum tensor from an isotropic state is determined by solving for AA*, AB*, and 
BB* , using the isotropic form at st * 0. 


“ 1^4 - kikj) 

Deissler (1970) has carried cut the analysis for the spectrum function 


E(k) ' k‘* 


^-(k/ko)2 


integrating (numerically) over wave space to determine the time history of the energy 
tensor. The solution depends only on the dimensionless time st and on the Reynolds 
number s/vk^'". For small Reynolds numbers, the energy decreases at all times, but 
for large enough Reynolds numbers the energy increases for a finite time and then 
decays. The y2 energy component decays at all times, as can be readily seen in the 
solution for v above. This in turn causes the cross correlation uv to decay, after 
some growth from its zero initial value, and the resulting loss of production (-8 Tiv) 
eventually leads to viscous decay. 


This linear picture of the development of shear turbulence is consistent with 



the experimental evidence at small times (st < 4) and predicts (surprisingly well) 
the magnitude of the shear stress correlation and the ordering 

> w"* > of the component energies. It has been determined experimentally, 

however, that the v component eventually grows, gaining energy through the 
pressure-strain correlation caused by the nonlinear tr^sport terms In the 

Navler-Stokes equations; this in turn leads to growth in *uv and the energy. The 

ultimate fate of the turbulence is currently a speculative issue, but there Is 
evidence of a self-similar evolution. The energy is still growing at the largest 
times observed In both the experiments and computation, and the ratios of component 
energies and Integral scales are still varying, although rather slowly. If the 
turbulence approaches a self-similar state with single velocity scale q and single 

length scale L, its attributes depend only upon t, q, L, s, and v. Dimensional 

analysis then requires 



If the velocity and length scales both grow monotonically at large 
times, sL/q must approach a finite nonzero constant, because if sL/q ->®, lir*^ ir 
theory is valid and predicts the ultimate decay of q; similarly, if sL/q-^0, the 
flow approaches Isotropy and q decays. With nonzero s the flow can not become 
isotropic in any case because for isotropic flow there is strong experimental 
evidence that L ' q^/e and q ^ t“^ from which it follows that d(sL/q)/d(st) > 0. 

Four shear flows were simulated, three of them (BSH9, BSHIO, and BSHLl) from the 
same initial conditions. The initial energy spectrum for these runs was simply a 
square pulse E(k) - 1 for 16 < k <32, a wave-number band containing roughly a 
quarter million Fourier modes. The initial spectrum for the fourth case (BSH12) was 
a square pulse at 10 < k < 20, the lower wave number allowing a higher Reynolds 
number to be attained. 

The shear rates and viscosity used are as follows: 


Case 

Spectrum 

(E(k)-l) 

s 

Viscosj 

BSH>» 

16<k<32 

20/2 

.OI //2 

BSHIO 

16<k<32 


.02//2 

BSHll 

16<k<32 

40/2 

. 02/ v'l 

BSH12 

10<k<20 

20 

.005 


The velocity-scale history of the four runs Is shown in figure 5(a). As 
observed experimentally, the energy decays at first but later grows monotonically for 
the remainder of the observed time; the growth beyond st-10 appears to be 
exponential. The g’uwth of the integral scales is shown in figure 5(b), with obvious 
anomalies beyond C'*10. The strange behavior, especially for BSH9, is caused by the 
fact that the two-point correlation is forming a rather broad negative region that 
diminishes the growth of its Integral because of the lengthening positive portion, as 
shown by figura 6. Similar behavior ic observed in the other correlations for run 
BSH9 in figures 7-9. By st-10, the largest scales of the simulation (rov».ghly the 
size of the computational period) have attained sufficient energy that they dominate 
the Integral scale. The evolution of the tht ee-dlmenslonal energy spectrum of run 
BSH9, which illustrates the point, is shown in figure 10. The sample in each 
spherical shell is proportional to k^ so that as the spectrum develops toward the 
large-scale end, more and more of the energy is contained in fewer and fewer Fourier 
modes. At st'lO the number of these energy-rich modes no longer provides an adequate 
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statistical sample I and although the simulation is still a periodic solution oi the 
Navier-Stokes equations, it is not representative, at the large scales, of 
turbulence* However the smaller scales still retain a good sample, and statiijtics 
from them (e.g., microscales) show no anomalous behavior* As shown in figure 11, the 
dissipation scales are adequately sampled throughout the run, but for the 

spectrum does not decay with wave number rapidly enough to allow complete resolution 
of the dissipation scales. All four cases appear to be adequately resolved at 
St 8; BSH9 and BSH12 are marginal at at ■ 10 \ and BSHIO and BSHll are marginil at 
st " 12. 

The primary dimensionless flow parameters are presented in figures 12 and 13* 
We have used the longitudinal integral scale in the streamwise direction (>: ) as 

the length scale and q * as the velocity scale* The mean shear s 

provides the time scale of the shear. In figure 12, the experimental data of TC 
indicate that the growth of velocity and length scales has equilibrated at st-S; this 
is prevented in the computation by the peculiar growth of the int€igral scale. The 
turbulence Reyr. Jlds number shown in figure 13 indicates that the experimental tlata 
are consistent with exponential growth In both velocity and Length scales, having 
roughly the same growth rate as the computation. The experimental Reynolds numbers 
are simply not attainable with the available computational resolution but we hope 
that the Reynolds numbers attained are high enough for the statistics of che 
energetic scales to achieve nearly asymptotic values. The classical mixing length 
(fig. 14) is formed from velocity scales of the turbulence and the (constant for each 
run) time scale of the shear; it appears to grow exponentially at j^arlier times than 
does q (fig. 5(a)). 

The distribution of energy among the velocity components is presented in figure 
15; it suggests the possibility that the component energies approach fixed ratios. 
If such a structural equillorium occurs between the tendency of the shear to produce 
anisotropy and the tendency of the turbulence to become ’.sotropic, the anisotropy 
measure of the equilibrium should depend on sL/q and (probably weakly) 
on qL/v. This Is shown more clearly in figure 16 where an adhoc linear scaling 

of the anisotropy is shown which collapses the data reasonably well. The discrepanoy 
between computation and experiment could be a Reynolds-number effect similar to tha: 
in the shear-stress correlation shown in figure 17. The Reynolds number shown there 
was chosen because it is the only dimensionless group independent of L and t. 
Although it seems physically reasonable to eliminate explicit dependence on time, the 
elimination of L simply allows us to avoid the use of the measured integral scales. 
The rapid rise of the shear-stress correlation from Its initial isotropic value of 
zero is evident. The shear -stress correlation is equally well collapsed by 
microscale Reynolds number. The reduction of the shear-stress correlation with 
increasing Reynolds number is primarily due to the fact that the contribution to uv 
of the small scales increases less rapidly with the increase of small-scale energy 
than does the contribution to the normalising factor (u^v^); this Is because of the 
tendency toward Isotropy of small scales. The explicit use of Reynolds-number to 
collapse correlations for medium and l.'irgc scales was first used by Stewart and 
Townsend (1951) in isotropic turbulence. In the shear flow, unlike isotropic flow, 
we have a direct energy measure, uv, that is relatively Reynolds number insensitive, 
and this can be used to advantage to normalize velocity autocorrelations. 

The autocorrelations presented in figure 18 are normalized by the single energy 
measure uv, and plotted against the separation, normalized by the single length scale 
the integral scale of Rn(r,0,0). The experimental points are from the 
Tavoularis and Corrsin data at X/H • 11. If the turbulence were in exact structural 
equilibrium this scaling would collapse all of the data, except near r - 0 where 
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viscous effects appear. At infinite Reynolds number the collapse would be complete. 
Perfect collapse re<}ulre8 that ratios of the component energies and ratios of the 
macroscales be universal constants. This Is too much to expect here^ however, for we 
have already st an that ratios of the component energies in both the experiment and 
computation are still varying slowly with sL/q and qL/v; moreover the 
integral scales of the computation are not as reliable as we would like because of 
the small sample of the largest energetic scales. 

The microscales measured by TC did not vary after st ■ S.65, but the computed 
microscales are growing slowly at st ■ 8, as shown in figure 19, with the growth rate 
diminishing with increasing Reynolds number. 


TABLE 1. - COMPARISON OF EXPERlMENTi^L AND COMPUTED TURBULENCE FIELDS 


Quantity 


u7^ /q^ 

u2^ /q^ 

up^/q2 


-uv/u'v' 


^11,2/^11,1 

Lii,3/Lii,i 

L22,i/I-11, 1 

L33, l/Lll , 1 
^12, ;/Li 1 , 1 
Li 2, a/i'll , 1 

^12/^11 

X13/A11 
X 21 /A 11 
X3 ]/Xii 

*^a/ *^c 
aa(‘) 


Vf/v * (-uv)/vs 
Rx * u'*ii/v 
sLii, i/q 
qLii.i/^ 

e/Lii.l = (-uv)/sLii,iV- 


Tavoularis- 

Case 

Case 


Corrsin 

BSH9 

BSH12 


X/H-11 

st*10 

at“10 

(Component energies) 

.535 

.508 

.485 


.186 

.185 

.209 


.279 

.307 

.306 

(Shear-stress 

.A5 

.49 

.49 

correlation) 

(Integral scale 

.33 

.44 

.59 

ratios) 

.25 

.17 

.14 

.23 

.41 

.33 


.34 

.23 

.17 


.90 

.88 

.88 


.40 

.69 

.74 

(Microscale 

.67 

.65 

.76 

ratios) 

.68 

.69 

.79 


.68 

.86 

.80 


.79 

.88 

.85 

(Principal stress 

4.3 

4.8 

4.0 

ratios) 

2.1 

1.9 

1.8 

(Principal axis 

20 

22 

24 

orientation) 

(Dimensionless 

179 

16 

38 

measures) 

266 

76 

104 


2.83 

2.23 

1.96 


3581 

219 

472 


.116 

.167 

.174 


A summary of the TC data at X/H - U (st'12.6) and of the computational results 
at st • 10 from cases BSH9 and BSH12 (our highest Reynolds numbers) is presented in 
table 1. The results of the computation agree well with the experimental data, 
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except for the Integral scale ratios (and quantities containing viscosity 
explicitly)* Part of this discrepancy is due to the fact that TC measured the scales 
by integrating to the first zero of the correlations and that we integrated over the 
negative portions as well; but this affects only l<n ,3 in the experimental data* The 
main cause of the discrepancy is that the integral scale » as stated by Batchelor 
(1953), is not as representative of the scales of motion containing the energy as we 
would like; it overemphasizes the largest scales where our sample is poor* In the 
experiment, a large sample of the largest eddies is swept past the probes by the mean 
flow, but in the computation we move with the mean flow and always observe the same 
large scales* 

There is some evidence in both the computational and experimental results that 
the turbulence approaches a structural equilibrium, with length and velocity scales 
growing exponentially. The experimental evidence, stated by Tavoularis and Corrsin 
themselves, although they seem to expect linear rather than exponential growth, is 
that the measured downstream transport Dq^/Dt and the production -suv of the 
turbulence are in a fixed ratio, within the experimental error, for 7*5 < X/H < 11. 
This implies that the dissipation rate is also in a fixed ratio to the transportT and 
this relationship was checked independently by measurement of the velocity derivative 
variances. This is equivalent to the fact that both the shear-stress correlations 
and microscales are constant. The energy growth rate is then proportional to the 
energy, and the growth is exponential. In table 2 we present a check of the 
experimental data for consistency with linear and exponential growth by using the 
data at X/H * 7.5 and X/H « 9.5 to predict the data at X/H * 11, Although either 
growth form extrapolated to X/H • 11 probably lies within the experimental error, the 
exponential form is consistently closer to the published values. The computations 
clearly indicate that the energy of periodic perturbations of a uniform shear grows 
exponentially, but the turbulence macroscale is bounded by the computational period. 
Unfortunately, by the time the computational velocity is clearly growing 
exponentially, the integral scales indicate an Insufficient sample at the largest 
scales, and the relevance of the exponential growth to turbulence (as opposed to 
periodic solutions of the Navler-Stokes equations) is suspect. 


TABLE 2. - DATA PREDICTED AT X/H * 11 
BASED ON DATA AT X/H - 7.5 AND X/H « 9.5 


Quantity 


Data predicted at X/H « 11 
Measured Linear growth Exponential growth 


u2 

.475 

.465 

.478 

vi 

.165 

.163 

.167 

w2 

.248 

.242 

.247 


.888 

.870 

.892 

I' n , 1 

57 

56 

57 


3.42 

3.26 

3.35 
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UNIFORM ROTATION 


The decay of turbulence in a region in uniform (solid-body) rotation presents a 
rather special problem to the turbulence modeler because in the Reynolds stress 
equations the net production term (production plus fast-pressure-strain) vanishes 
when the turbulence is isotropic. This is an artificial situation in the sense that 
it is probably impossible physically to generate isotropic turbulence in a fluid 
rotating as a whole, but mathematically it is simply a matter of specifying isotropic 
Initial conditions. At the Reynolds-stress level it is not possible to determine 
whether initially isotropic turbulence will remain isotropic. The experimental and 
computational results indicate that the turbulence does become slightly anisotropic, 
and this is consistent with the linear theory valid low Rossby number. 

At small Rossby numbers (q/i^L) the mean rotation (Q) causes plane waves 
e^!S*5 to propagate with phase speed 2fi*k/k, giving the Fourier modes 
some physical ^*wave** significance (unlike most turbulent fields in which the wave 
space is purely a mathematical convenience). The nonlinear t^.rms can then be viewed 
as direct wave interactions. The bilinear convection term in the Navier— Stokes 
equation is interpreted in wave space as the interaction between wave vectors k' and 
k” to give the convective contribution to the time derivative at wave vectoF k^+k*’ 
so that the time derivative at Ic is determined by the convolution sum over all wave 
vectors k", k" satisfying k"+k*‘ ■ k. Taking the y-axis as the axis of rotation, the 
Fourier modes at low Rossby number have the form 


rUk.t) - A(k. ^ ^),-l(2k,/k)St 


This is a classical ”two-tlme-scale" problem at low Rossby numbers; the 
"long-time” is qt/L where q and L are representative of the turbulence and the 
"short-time” is fit. We are interested in the long-time variation of the spectrum 
tensor, and this requires that we integrate the wave interactions over cimes (T) 
much larger than the short-time scale but much smaller than the long-time scale. The 
correlation over this intermediate time between u(k) and the contribution to its 
time-derJvatlve due to the Interaction between u(k') and u(k") is proportional to 


sin u)T 
uT 


where 


00 




For very small Rossby numbers, only interactions between waves k', k" having 
kj/k ± kf/k' ± k' 2 /k" « 1 are significant on the long-time scale. Rotation then 
acts to attenuate the effectiveness of the nonlinear interactions, which in isotropic 
flow move energy to higher wave numbers and enhance the dissipation; this explains 
the experimentally observed reduction (Wigeland and Nagib (1978)) of the dissipation 
from its value without rotation. This argument only requires that rotation attenuate 
interactions and assumes nothing about which directions in wave space or which 
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velocity components suffer the reduced energy flow. However the experimental data of 
Wigeland and Nagib (1978) indicate that the energy is redistributed among the 
velocity components with the ratio of energies reaching an equilibrium when the 
transverse macroscale begins to grow at a faster rate than the longitudinal 
macroscale. 

We have simulated the flow at five rotation rates, but with a microscale 
Reynolds number of only about 3, because of a program input error. As a result, 
dissipation dominates the Reynolds-stress equations. • The dissipation rate la 
nevertheless reduced slightly (•■25 percent at the highest rotation rate) at early 
times, as shown in figure 20. A structure parameter indicating component energy 
anisotropy Is shown In figure 21; the longitudinal component W gains energy through 
the slow-pressure-strain term. As the tables in appendix A indicate, the 
contribution of the fast-pressure-strain term varies widely in magnitude and takes 
both signs, Indicating that it Is oscillating. The slow-pressure term, on the other 
hand, varies slowly and alwa‘/s works to enhance The experimental data show 
higher anisotropy, but this could well be a result of the anisotropic initial 
conditions of the experiment or to the Reynolds number disparity. 

In figure 22 the anisotropy at the dissipation scale is also shown to reach an 
equilibrium, and the correlation with fit, rather than vt, indicates that this is 
not entirely a viscous effect. 


TURBULENCE MODELS 


The results of the numerical simulations can be used to aid the construction of 
turbulence models at any level, from simple algebraic models of the Reynolds stresses 
to spectral models such as the quasi-normal approximation. Appendix A contains 
information at the Reynolds-stress level for a sample of the computed fields 
discussed in the previous sections. In this section we illustrate how these data 
might be used to test proposed models or to determine "model constants" for some of 
the modeled terms in the Reynolds-stress equations. 


THE REYNOLDS-STRESS EQUATIONS 


When the turbulence field is homogeneous, the Reynolds-stress equations become 

dt “ "Ui.k^juk - Uj,k<Jiuk + 2psij - 2vui,kuj,k 

, I pressure- .. 

production . disslnatlon 

strain 


where the mean flow gradient Uj[^^ dependr, only on time, and is the turbulent 
strain-rate tensor Sij » j+uj ,i.)/2. The over-bar denotes an average over physical 
space. Pressure Is decomposed into its so-called "fast" and "slow" parts given by 
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t iaite: [iiiflUliHf ' 




(slow pressure) 


7?p(l) 




V^p(2) - -2UijUj j 


(fast px'essure) 


The terms "fast" and "slow" refer to the fact that p(2) depends directly on mean flow 
gradients and can therefore have Its time scale imposed directly by the mean strain 
late.^ The slow term, on the other hand, responds indirectly through changes in 
turbulence structure. For homogeneous flows we have exactly 


2p 


ij 


a'Pf U 


e,m 


where 


2 2 --y «i(b)^m(b) 

The sum indicated is over all Fourier modes. The fourth rank tensor a has the 
obvious symmetry “ 


oiiil _ ,im _ „mi 
ej ‘^.1 


and the continuity condition, k|U£ * 0, implies 


mi _ 

ei 


In addition, the contraction ajj - 2uiun, provides a convenient scaling factor, 
and as a factor of the mean deformation rate it emphasizes the close relation between 
the production and fast-pressure-strain terms in the stress equations. In essence, 
the production and fast-pressure terms sum to produce the "net" production. Tables 
of the tensor a, shortened by use of Its symmetry, are Included for each recorded 
field in appendix A, with the budget of the stress equations above. 

In the Reynolds-stress equations only the production term can be computed 
directly from the stress tensor. The dissipation and fast-pressure terms are 
second-order velocity moments like the stress; however, thev contain velocity 
derivatives and thus depend on the spectral distribution of stress (the spectrum 
tensor) and not merely on its average value. They must be modeled along with the 
8 low— pressure term, which is a third-order velocity moment. 

We consider first the slow-pressure-straln and dissipation terms, following the 
analysis of Luraley and Newman (1977) who combine the slow-pressure-straln tensor with 
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the deviator of the diesipation tensor* Lumley writes 


2p 


( 0 , 


ij - f 3 




where ^ is a traceless tensor and e is the dissipation rate of kinetic energy 




which must be modeled. The approximation (Rotta (1951)) that ^ depends only on the 
tensor b, 


K 1 X 

- 3 ^Ij 


and on scalar attributes of the flow, requires for proper tensor invariance 
that ^ be an isotropic function of b, 


♦ ij “ ‘^‘>ij ^(bikbkj + ^ 


The coefficients B and y are functions of scalar invariants of the field such as 
Reynolds number, q^/ve, II, lit, 


where -2II * b^jbji , 3111 » b^jbji^b|^i 

are the Invariants of b (the third independent invariant, is zero by the 

definition of b). Lumley and Newman show that for high Reynolds numbers , Y“^» 

slightly anisotropic flows are to return to isotropy* 

We will retain for our purposes here only the linear term of the model* Some 
information about the coefficient of the linear terra B can be determined by 

considering limiting cases. In the absence of a mean flow (for example the "return 
to isotropy" following an applied strain), B 2 If the anisotropy is to decay* For 
any anisotropy as the Reynolds number vanishes since in that case the velocity 

components all decay in proportion to their magnitudes and the anisotropy b is 
constant. Lumley and Newman fit the slightly anisotropic data of Comte-Bellot and 
Corrsin (1966) and find that, again, B *^2 as the anisotropy vanishes at high Reynolds 
number in the absence of a mean flow* In later papers, Lumley (1978, 1979) shows 
that In order to force the energy tensor predicted by the model to be realizable 
(have positive definite component energies in principle axes) S^2 as a component 
energy vanishes. Lumley found a function of the invariants 
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■S:*r.n.i;Trr = 








iin 




G(11,I1I) - 1/9 + 3111 + II 

that vaniahea If and only If a component energy vanishes , and he suggests that this 
function be explicitly Included In the model to ensure realizability. For example, 

B - 2 + B'fidl.III) 


where B ' is a function of scalar invariants. 


The relaxation toward isotropy of turbulence subjected to axlsyrametric strain 
has been simulated, taking as initial conditions fields at total strains 2 , 2, and 

H from the set of axlsymmetric strain runs discussed earlier. In figure 23 we show 
the variation of ^ and its pressure-strain and dlasipatlon terms with the anisotropy 
b. In this flow only the normal stresses are significant, and two of them are nearly 
equal, the deviation being a result of the finite sample. The two families of points 
in the pressure-strain plot correspond to the two values of viscosity in the 
simulations, indicating the Reynolds-number dependence of the coefficient B. This 
Reynolds-nuraber effect is more apparent in the dissipation and pressure-strain parts 
than in their sum indicating the advantage of treating these terms together. 

This is done by default when experiv.\entai data are used to determine model constants, 
because the dissipation is assumed to be isotropic in order to find the 

pressure-strain (plus dissipation devlatur) from the Reynolds-stress equations. The 
model coefficient 6 for each field was found by least-square fitting the 

measured ^ to the measured The fit is nearly perfect for each field. Inis is 

expected since i and b are both diagonal, traceless, and axlsymmetric and thus can 
be .'elated exactly by ”a single coefficient. The resulting coefficients for all of 
the fields are plotted against a Reynolds number q^/9ve and G(II,T1I) in figures 
and 25. Although the scatter in the coefficients is large, the trends are consistent 
with Luinley's predictions. The Rotta model applied to the axlsymmetric strain runs, 
during the straining period, gives essentially the same results. 

In figure 26 we present the data from the plane-strain cases discussed earlier, 
the groups of four points being associated with the four strain rates applied. The 
single degree of freedom of the linear model does not fit the two degrees of freedom 
of the data well, but the fit Improves as the strain rate decreases. The nonlinear 
model, having two degrees of freedom, would fit the data exactly. 


The anisotropy of the shear cases Is plotted in figure 27. The fit achieved by 
the Rotta model is shown in figure 28 as a plot of the measured ^ for each field 
against the predicted for that field by the linear model, using the measured b 
and the coefficient B giving the least-square fit for that field. Although this 
would appear to be a more difficult case than, the two-dimensional strain, since four 
elements of the stress tensor are active, the model seems to fit the data fairly 
well. The coefficient grows with Reynolds number (fig* 29) as predicted by bum ey, 
but its variation with C(TT,m) is uncc-tai i (fig. 30). 


The modeled terms for the unifoim rotation cases are presented in figure 31. 
The dlssipc-itlon anisotropy is correlat! ,1 with the stress anisotropy, but the pressure 
strain is negatively correlated with it. The pressure-strain term Is evidently the 
source of the anisotropy in this case, contrary to Its usual role as the "return to 
Isotropy" mechanism* 

The performance of the Rotta model must be judged by how well it 
elements of the tensor ^ at each field and how well the coefficient for each field 
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can be fitted Co the scalar invariants of the field. With the exception of the 
plane-strain runs, which seem to be sensitive to strain rate, the linear model 
provides a fairly good fit to the data with coefficients Chat vary (qualitatively) 
with Reynolds number and anisotropy level in accordance with Lumley's predictions. 
The results do suggest, as has Lumley, that the time-scale ratio between the 
turbulence and mean deformation be included in the list of invariants determining the 
Rotta coefficient. 

It would seem that modeling the Censor 

- 2 E ^ 

iiC 


would be an easier task than modeling the slow-pressure term which is a third-order 
velocity moment; however, this has not proved to be the case. The most widely used 
model of the tensor a is that of Launder, Reece, and Rodl (1975) (LRR hereafter) 
which has the form 


a 


mi 



+ CB 


mi 


where C is a coefficient to be determined (again, a function of Invariants) and A 
and B are fourth-rank tensors of known form that depend linearly on the 
Reynolds-s tress tensor, satisfy the symmetry required of a, and vanish when 
contracted over ij . The contraction over Ij gives 


. 23pr , B"! . 0 


When the turbulence Is isotropic 




e i 


. 2E 






in which case 


ami 


0 


and 


mi IB ■ 

'ij ''ij 
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mi ij 


mi ij 


^mj^ii 


) 


which gives the exact result for an isotropic flow regardless of the value cl C • 
The model has certain fallings, as do all models. Some of these have been pointed 



out by Lumley (1978), who shows that the model can not meet the reallsabillty 
condition, and by Leslie (1980), who demonstrates that the model is too simple in 
form to represent the general fourth-rank tensors a that might occur. This latter 
argument is easily verified by noting that regardless of the coefficient C , the 
model produces the following spurious symmetries (no summation implied) 


a 


a 


a 


ki 




ii o 


ij 

ji 

- 


a 

ii 


i ^ j 7^ k 


i j 


The addition of nonlinear terms (higher-order products of the tensor) would remove 
these spurious symmetries and introduce more coefficients, and as noted by Lumley, 
might allow the model to meet the realizability constraint* The .justification for 
nonlinear dependence on the stress-tensor, when a is linear (formally) in the 
spectrum tensor, is that a depends also on the distribution of energy over wave 
space, a dependence lost when integrating the spectrum tensor to the stress tensor. 
The idea that the anisotropy of the spectrum tensor with respect to k might be 
expressible in terms of the anisotropy of the stress itself leads to the possibility 
of nonlinear dependence on the stress tensor, and, as suggested by Lumley, the 
connection might be made using the theory of rapid-distortion. In view of the 
difficulties associated with the LRR model, several authors (e.g. Leslie (1980)) have 
suggested modeling the rapid-pressure-strain term, or even the total pressure-strain, 
as a unit. This avoids the detailed constraints on the model imposed by the formal 
definition of a, but the ability tr handle flows with large variations of imposed 
mean strain would be lost. 

We compare the results of the LRR model with the computational data for the 
shear runs in figure 32. The 36 different elements of the tensor a are each 
normalized by q^. The constant for each field, found by least-square fit to the 
data, is plotted against Reynolds number in figure 33, and against G(I1,II) in figure 
34. The model succeeds rather well in fitting the 36 different elements with a 
single constant for each field, and the constant appears to be related to the 
anisotropy level, If not to the Reynolds number. 

The least-square procedure, which was used simply to Illustrate how well the 
models fit the data, indicates that there is room for improvement. The use of the 
data to determine the coefficients of a postulated model can be done in a number of 
ways. For example, rather than finding coefficients for the Rotta and LRR models 
separately, they can be found simultaneously by fitting the sum of the individually 
modeled terms. This is how experimental data are used for the shear flow (Leslie 
(1980)). When this procedure is used on the computed results, model coefficients 
very close to the “recommended values" are found. 

The real value of the data of appendix A howv^ver, is not that it allows model 
coefficients to be determined in a more precise way than does experimental data, but 
that it provided the detail needed to determine the range of validity of a model and, 
it is hoped, to suggest how the model might be modified to increase that range. 


CONCLUSIONS 


The large-scale direct aimulatlon of homogeneous turbulence begun by Orssag and 
Patterson a decade ago, has advanced wlrh increases in computer power to the point 
where today it is possible to capture the statistics of the energetic scales of 
motion at low Reynolds number for mouerate imposed mean deformations. Further 
advances in hardware and method will permit higher resolution in the future, but the 
resolution available today is adequate to allow the extraction of information useful 
to both model builders and **pure” theoreticians* 

The relevance of the aimulatlon to real turbulence has been established, in this 
author s opinion, to the point that the results can be used to fill in the gaps in 
the experimental database, and It Is hoped that the information Included for this 

purpose in appendix A will lead to closer scrutiny of the fundamental assumptions 
implicit in turbulence models* 


t 
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APPENDIX A 


LISTING OF DATA FROM A VARIETY OF COMPUTED TURBULENCE FIELDS 


The *‘raw data‘* of a variety of computed turbulence fields is presented here in 
tabular form# Fields from the following (deforming mean) runs are included* 


Deformation Rate 

(r^T 

(Plane strain) 4.5 

-(dv/dy)l ® 

d«/dd ) 


-(dv/dy)l 
dw/dz i 


(Axisymmetrlc 

contraction) 


du/dx 
-2 (dv/ dy) 
-2(dw/d?0 


BSH9 

BSHIO 

BSHll 

BSH12 


(Shear) 
du/ dy 


(Rotation) 


du/dz 
- (dw/dx) 


Viscosit) 

■02f=Tf 

. 01//2 

. 01//2 

.OI//2 

.OI //2 

.OI //2 

.Olh^ 

.O 2//2 

.O 2//2 

. 02/^"2 

.OI //2 

. 02/*'2 

.O 2//2 

.005 

. 01/»'2 

.Qlh^ 

.01 

.Olh'l 


Inltic^;! Condition 

(Developed 

Isotropic 

field) 


(Developed 

isotropic 

field) 

(Developed 

isotropic 

field) 


(Undeveloped 

isotropic field, 
square spectrum) 


(Developed 

isotropic 

field) 


Fields from these runs are named using the run name above followed by a letter 
designating the individual field. The plane strain fields are given at nominal 
strains of 2, 2/2, and 4, with suffixes A, B, C, and D, respectively. The 
axisymmetric-strain fields are given for nominal strains of ^/I, 2, and 4, with 
suffixes A, C, and F, respectively. For the shear runs, the suffixes are not 
consistent from run to run, but for all runs the fields are given for st “ 2, 4, 6, 
8... The rotation fields are taken more or less evenly in time and are labeled 
sequentially A, B, C, D, E, and F for each run. 

Runs named RX... are relaxation-to-isotropy cases for each of which two fields 
are given; for example, RX24A2 and RX24A4 are fields from the run using as initial 
conditions the turbulence field CD24A but setting the strain rate to zero. Thus the 
data of field CD24A with net production set to zero provide a third field of the 
relaxation run. 




The following data are presented for each field: 

T ■ time with units [T] 

Vise ■ V (Kinematic, viscosity) 

RII ■ ujui = (Velocity variance • [L^T“^] 

2*kinotlc energy) 

PlI ■ 2U^ jUji^Uj (Net production of RII) [L^T 

Dll •• 2e « 2vu^ j (Net dissipation of RII) [L^T"^] 
IX ■ bijbji (Invariant of anisotropy tensor) 

III “ (Invariant of anisotropy tensor) 

MEAN DEFORMATION DU/DX etc. [T“M 


INTEGRAL SCALES 



(Note the absence of 2's in these length definitions.) 


REYNOLDS-STRESS BUDGET 


IJ - ij 


Stress component 

RIJ/RII 

- ^i^j/^k 

Component energy ratio 

RATE - 

Id, . 

27 ^ 

Net time rate 

PROD - 


Production 

P2S - 


Fast-pressure -St rain 

PIS ■ 


Slow-pressure-n train 

DISS « 


Dissipation 
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...^402 9.0411 

,0174 -.3004 

0.0001 -«om 

8.C0U 0.0011 
• .C004 -.4002 

C.COCO 4.0044 


22 

0.1174 
-.6001 
O.OCOI 
0.1114 
0.001* 
0 . 40 40 


0.0004 

-.1*04 

-.0014 


. 01V9 

-0.17* 
-0.001 
-0.000 
•4.411 
0.001 
-0.413 


n 

0 . 1*40 
--000T 
-.0010 
0.914? 
0.0444 
0.1144 




Vise 

Rtl 

Ml 

OH 

It 

III 


1.204C-02 
T,0TM-01 
9.993C 00 
4.1221 01 
4.1141 01 
2.2101-02 
-1.394 1-01 

IMTtORH SCALIS 
X Y t 

*1001 .0*13 .091* 

,1021 *IC4« .02*4 

.1009 .6277 .192* 


RCAh OIRORRATICR RRTf 
I Y 2 

40.00 0.00 0.04 

o.eo -20.00 0.00 

0.00 0.00 -20.00 


,0721 .0132 *0422 

,«*Tl .0*17 .0341 

.CkT) .0)42 .0414 


Rt«/Rll 
0.21 1 
•0.000 
0.941 
0.341 
•0.000 
0.141 


RCYTiOlOVItOM) lUDOtT 


RAM 

•0.414 

0.000 

-0.001 

0.243 

-0.(«fil 

0.242 


RROO 
-1.011 
0.040 
-0.001 
0.49 7 
-O.OOl 
0.1*9 


*79 

0.91* 

•0.001 

0.001 

•9*219 

o.oco 

-0.272 


*19 

0.491 

0*401 

•0.000 

-0.01* 

9.000 

-0.0)9 


tJ/RI 
It 0.071) 
12 -.0000 
0.0001 
.17*9 

9*0004 
.1791 


*Rl95yRI-StRA|R TIMOR 


IS 

22 

21 


12 

--0C(2 

•*ni*o 

0*0004 
O.OOOl 
-.00(1 
- .OOCI 


11 

C.0007 

-.009} 

-.0972 

0.0004 

-.0000 

0.0004 


22 

0.1)44 

-*0001 

G.O804 

0.1712 

0.0400 

0.41*1 


13 
*.0002 
-.9000 
0.390) 
O.OMO 
-.1392 
-.0004 


> *13S 

-0.114 
0.000 
'O.OOl 
-0.319 
'0.001 
-0.3J* 


33 

0.1744 

-.0013 

-.0002 

0.411* 

0.0003 

0.1724 


t 2.4)97-02 

Vise 7.0714-0) 

RII *.1111 40 

*11 1.4771 02 

mt i.;ia 02 

n 4,43)1-02 

til -1.29*9-02 


RUN 0C7CARAT1CR RAM 


40.80 

0.00 

4.00 


0,00 0.00 
-20.00 0.00 
0.00 -20.04 


IHTIORAW UAlM 
I Y 2 

.0*02 .0401 .09*4 

.2011 .0113 .4102 

,203* .0U3 .01)* 


FtlCRCSCAklS 
X V 2 

.0444 *0)** .03*4 

.1294 .9922 .0)04 

,12*4 .9)10 .0924 


RtJ/Rli 

0.077 

-0.000 

0.001 

0.4*C 

-0.001 

0.4ft) 


*|YRU 09-9TAi:S twoott 
AA1I • *400 * *29 • 

..3.249 6.12* 

0.000 
•O.OOl 


-0.1*9 

'O.OOC 

•c.eoo 
0.122 
•0.002 
0. 124 


o.irT 


-0.001 

O.OCl 

-0.04) 

•0.000 

-0.0*1 


*IS 

0.0*4 

.'•000 

.>.000 

0,0)- 

O.OO- 

-4.01 


.J/RI .. 
11 0.0221 
12 -.0901 

1 ) 0.0001 

72 0.0«»7 

.3 0.0002 

33 o.e*«4 


*FAST> PRE1SU4I-S1R41* MNSCR 


12 

-.0CC2 

-.0111 

O.OOOl 

-.0000 

-.00C9 

..4:(7 


13 

0.0001 

-.0002 

-.0119 

6.0001 

0.0001 

C.0804 


22 

0.0)9* 
o.coco 
Q.OOOV 
0.21*1 
0.CO13 
Cl. *299 


2) 

0.0904 
• .0002 
4.3002 
-.0009 
-.2071 
-.0014 


. OtSS 
-0.049 
0.000 
•0.000 
-0.494 
0.401 
■8,497 


9.0091 
C.*)20 
0. 9*)4 7 
0.2U't 


• •(•4v»*»«* 44M4 — fO)99 

T 1..'l92-0i ' 

Vise T,6.’lf-0) 

Rlt 1.4*. t 30 

*n 3.41)k 00 

01) 2*0141 41 

1 1 1.4701-1.7 

III -1.121I-01 

|tlM64Ak scuts 
» Y I 

U .111) .0*79 *0*T» 

V .1192 .1497 .6423 

* .1047 .0441 .1400 


■ t Ih 0CFCR4.4) ICR R4*£ 
A r I 

9.00 9,00 0.00 

0.00 -2.9M 0.3- 

0.00 9.00 -7.90 


.0422 .0931 .0132 

,(794 .041* .0924 

.0741 .09)0 ‘OMt 


AIJ/Rli 
0.214 
0.000 
0.901 
0.)4 7 
-0.001 
e. )4i 


4(<iFiCl.0S SIRISS IV06I7 


RATI 

-0. )(7 
O.OOl 
•O.OOC 

-C.2I0 

o.coi 

•0.2(9 


*4.VS 
-0.2.9 
-0.004 
-0.000 
0*141 
-0.091 
0. 144 


*29 

0.127 

•0.000 

0.001 

-0.0*2 

O.OOC 

-.',0*9 


*IS « 
9.042 
O.OOl 
• O.OQ*] 
-0.020 
O.OOl 
•0.022 


i2/A| 11 

It 0.04*2 

W -.0000 

1) 0.0004 

72 0.1717 

2) 3. 0010 
t> 0.171* 


•fast* FRtSS'jll-STRA]. TIMOR 


12 

-.0601 
-.042 ! 
0.0409 


1) 

v.ecc* 

-.3002 

-.0440 

0.0412 

-.ficei 

O.OOC) 


0.21*4 

-.0002 

0.0004 

0.1A17 

0.0004 

0.))4* 


13 

-.000) 

-.oo'?:' 

.'.3004 

-.1007 

-.K49 

-.00.) 


0199 

‘0.2*1 

0.000 

•0.000 

-0.94* 

0.000 

•0.)7l 


)1 

0.2124 

-.0040 

-.004* 

0.4017 

0.0004 

0.1*29 


fill 

HI 

*11 

Oil 


2. ))2l-3t 
1.0719-03 
l.Q*lt 30 

3. *371 00 
1.*04l 00 
T.7A9C-02 

-*.•)« 1-4} 

IMIORAl 1(A179 
I 1 ) 

.0010 .1007 .04** 

.24*4 .17*2 .3114 

.237* .012* .142<> 


*IAIt CMC4RATKR *At2 
I Y 2 

1.00 3.00 O.OO 

0.00 -2.90 O.OO 

3.00 0.00 -2.90 


.1114 .Ofttd .041* 
.1970 .1044 .0*94 

.117* .0*97 . 1047 


*!YRCL01-J)«M> )W0C»» 


9.19* 

0.909 

3.C31 

0.445 

-0.070 

0.4*1 


44tl 

- 0.121 
•0.001 
-0.90) 
-C.ll > 
O.OOC 
• O.UJ 


*400 

-0.201 

-o.oor 

'0.001 

0.4)1 

-0.000 

0.42* 


*2) 

0.121 

‘0.041 

O.C«l 

-Q.C»2 

•0.001 

O.Ot) 


*1S 4 

6. oil 
0.001 
•0.60* 
-0.044 
0.902 
-0.044 


U2«l U 

11 0.04)1 

12 '.0904 

\) 3.0001 
22 0.0131 

2) -.0900 

13 9.0IM 


•*ASt* F4fSSK«l 9>»41» TIMOR 


U 

• 30C » 

• .0214 
O.OCOI 
o.ccci 
- .00(1 


23 
0.0004 
0.0011 
O.OOCI 
0.0004 
-.17*4 
-.002) 


3193 
•0.14) 
-O-OOl 
■0.300 
•.*.42 4 
O.ObO 
-o,4;« 


1) 

1.1401 

0.<)00l 

-.0911 

0.9900 

0.0001 

0.1140 


»*4t*4«***»MM1 


1.4932-02 
7.on*-os 
1.4141 00 
1.293{ 62 
1.002t 02 
9.1421-02 
-6.I14»-63 


Vise 

HI) 

*lt 

on 


t V***************** 

ream DEEORNATtCR RAM 
X T 2 

40.00 C.8<i O.OO 

0.00 -20.00 6.90 

O.VO 0.00 -20.00 


INTCORAl SCALES 
X V 2 

.0441 .0934 .09)0 

.13*4 .1(14 .0144 

.1144 .0141 .04)7 


RICROSCALIS 
X Y I 
,0741 .040* .040) 

.94*0 .0974 .9392 

.044) .4391 .9140 


R1J/4I1 
0. 147 
-0.000 
0*001 
0.42 9 
-0.000 
0.42 4 


riynolos-itress iuoget 


RATI 

-0.4*9 

-9.000 

-0.001 

0.34A 

-9.002 

0.344 


RROO 
-0.4*1 
a.MO 
•0.001 
).*92 
• >.001 
1.494 


*29 * *19 
0.321 0.094 

-0.901 0.000 
8,001 0.000 
-O.lll -0.024 
0.000 -0.000 
-0.14! -0.024 


LJ/NI 11 

11 0.04TT 

12 -.0040 

13 0.0001 

22 0.1233 

23 0.0009 

)) 0.1234 


•PAST* *R|SSURf-STR*IR TIMOR 


12 

-.0092 

-.023* 

0.0603 

0.0001 

-.0006 

-.OCCl 


O.OOOl 

-.000) 

-.0141 

0.6001 


0.133T 0.0001 

-.0601 -.0001 
0.0004 0.0993 
0.1424 -.OQCl 
-.1*44 
-.0004 


> 0IS9 
-0.1*3 
0.009 
-O.OOl 
-0.414 
-0.001 
•0.414 


32 

0.1242 

-.0019 

0.0091 

0.9249 

0.9004 

0.14)0 


3.91TE-02 
,C 7.0TM-02 
A.4541 00 
2.1 T2E 02 
1.4*4C 02 
I.M9T OJ 
I -1. *221-02 


t ••••••«••••••*•••• 

RIAN OtRCRRATKR RATI 
ATI 
40.00 0.00 0*00 

0.00 -20.00 0.04 

0.00 0.00 - 20.00 


INMGRU 9CAlC:S 
X T I 

.0T2* .0(42 .094* 

.29*4 .<19* .0077 

.2902 .6(79 .0402 


RICR09CRLE9 


, 1006 .0)44 .0314 

.192* .0102 .0249 

.1912 .0244 .0904 


RlJ/RIi 

0.093 

•0.000 

O.OC-'. 

0.47. 

•0.001 

0.47) 


REYNCLOS-StRISS lyOOIT 


RAM 

-0.117 

-0,000 

-O.iSOO 

0.211 

-0.002 

0.304 


*ROD 

-..,114 

0.000 

-O.OOl 

0.434 

-0.002 

0.04) 


*29 

0.07* 

-0.001 

0.000 

-O.Cll 

-0.000 

-0.034 


*15 

0.0*2 

0.000 

0.000 

-0.030 

0.001 

•0.931 


LJ/NI 
11 0.014) 

i; -.0001 

1) 9.0001 

22 0.047* 

2 ) 0.0001 

1) 0.047* 


-FASf *RMSVRI-STRA|R TIMOR 


IS 

-.0062 

-.QC72 

0.0000 

-.00(1 

-.00C4 

-.0062 


1) 

O.OOOl 

-.0002 

-.0071 

0.0007 

0.0001 

0 . 000 ) 


22 

0.0940 

O.MOO 

0.0604 

O.IaM 

0.0017 

0.**0) 


23 
0.0009 
-.9001 
0.0002 
-.0007 
-.1141 
-.0017 


* OtSS 
'0.040 
•0.006 
- 0.000 
•0.470 
•O.OOl 
•0.470 


33 

0.09*2 

-.9012 

0.0001 

0.4*40 

0.0004 

0.2270 




1. 3411-01 
r.etiHO) 
1.44*1 60 
3.«0)t 00 
1.1226 01 
4. 1214-62 
-3.4146-0) 


Vise 

Rtl 

*11 

on 


RIAN OIRi»RA''ICR RAM 
* T 2 

) .CO 0.48 0.00 

6.00 -2.90 0.00 

0.00 0.00 -2.90 


INMOAAh SCALIS 
X V I 

.101* .0111 .071* 

.191* .1*39 .0)19 

.14*4 .0128 .19(1 


RieXOSCAltS 
X V 2 

.0410 .0974 .0940 

.0443 .0430 .6974 

.0417 .8974 .0411 


RlJ/RII 

0.1*7 

0.091 

6.091 

9.414 
-O.t 11 
0...14 


RIYNOL09-STXEM BuOOII 


RAM 1 

-0.239 

- 0.000 

-o.ooo 

-0.214 

6.002 

-0.22» 


RROO 

-8.214 

•0.000 

- 0.000 

0.2*7 

* 0.001 

0.271 


*2$ 
ft. 134 
•0.001 
0.001 
-0.C46 
• 0.000 
•0.0*4 


PIS 

0.04* 

0.061 

• 0.001 

'0.031 

9.002 

•0.0>* 


LJ/Nl 

11 0.0*42 

12 -.9Q«^ 

13 8.0003 

22 0.1)29 
21 0.000) 

23 0.1324 


FASf 


12 


.'.(>006 

-.6)^*' 

0. 0QCl 

1. CU0 
-.040.: 
O.vCCS 


;.06U 


PR tSSWRi - strain MRSO.1 

22 2) 
0.1991 0.0900 

-.0063 
e.cccft 
0.1791 
0.0012 
0.49*2 


-9)92 

o.ooi: 

-.0901 

C.0092 


6.9002 

9.000) 

-.0009 

-.1414 

-.0019 


QISS 

-6,217 

‘•9.090 

- 6.000 

-0.114 

O.OOl 

-6.341 


19 

9.199* 
-.0011 
-.0014 
0.4*92 
0.OC04 
0.17 7'' 


T I.TIlt-01 

Vise 7. cm 6) 

Rlt 4.44*1-01 

*11 3.*4*f 06 

on *.7)ftl 90 

n 4.1101-02 

III -1. 12*6-02 


RUN enORNAMCM RATI 
A T 1 

9.00 0.06 6.00 

0.00 '2.36 0.00 

6.00 0*00 ‘1.90 


INTICRAL scalis 
I Y 1 

.0/** .1CI4 .1004 

.3014 .1)42 .0197 

.24)4 .0274 .1444 


HtC40SCALES 
X Y f 

.193* .0711 .0M4 

.1419 .1134 .0*72 

.1424 .0*74 .1143 


RU/RIl 

t.otr 

t.OOS 

O.OOl 

0.499 

- 8.001 


*t7NOL09-SfR|SS Iu06lt 


RAtl - RROO 
-O.C44 -0.1*3 

0.006 -6.000 
-O.CyOl -9.000 
-0.041 0.4*0 

-0-862 -4.081 


8.49* -0.9*1 9.1A3 


■ *2$ 
9.107 
-9.001 
0,001 
-8.613 
-8.402 
-0.649 


. *19 

0.949 
0.002 
-9.801 
•8.844 
9.801 
‘0.04* 


LJ/Nt 

0.0)34 
-.0804 
8.8002 
8.6*47 
0.6900 
0.8447 


OFASt* *RI9tuRI-SM41R UN90R 


It 

1) 

22 


12 

-.40(2 
-. 01 ** 
-.0681 
8.60* 
-.0689 
- .000) 


13 

0.0099 

-.0604 

-.8171 

9.661) 

-.4602 

6.9440 


12 

0.1124 
0.9001 
6. MCI 
0.26*4 
6.0024 
9.9404 


I) 

0.9010 

0.691* 

-.4001 

9.6094 

-.1406 

-.0921 


• OISl 
•0.117 
- 0.001 
- 0.000 
-0.412 
- 0.001 
-0.441 


3) 

0*11)4 

0.0094 

-.6014 

0.9494 

0.9049 

0.207) 


25 




r 

viK 

mi ».7m oe 
mi ».2iof 
on 4.1111 01 

I I *.2241-0) 

III -l.Om-04 


fl«K Stromi4YICN 04tf 
( T < 

10.09 O.OO 9>09 

0.00 »**w o**"* 

0.09 0.S9 >9.09 


INtfOMAl )C4LL'S 

I y 2 
,!)«/. .0i«<l .0704 

.0414 .1444 .0912 

.0441 .OUl .11*2 


.1004 .04)1 •04)1 

.0024 .9414 .0941 

.00)2 .0944 .9424 


rnrhCkOi^ifoios ouoort 


4000 • 02$ 

-0.941 C.141 

>0.900 0.009 

-0.000 0.001 

0.2)0 -0.(44 

-C.OOO 0.090 


02$ * MS < 

C.141 0.004 

0.000 0.000 

0.001 0.090 

.0.(44 -0.002 


’f4$I» OOf >SUO«>)tO*J* flhSO* 

LJ/NI 11 12 1) 22 2) 

11 9>I00(. O.COOO O.COIO 0.22)0 -.0004 

12 9.09C2 -.04$$ -.0004 -.900) 9it90»9 

1) 0.0001 0.90C7 -.9)11 *.90(9 »‘9«9» 

22 9.2114 0.<j0e9 0.9099 9.1962 0.9001 

2) 9.0001 -.9(14 


909) 9^9099 

90(9 9.909) 

1962 9.9991 


)1 0.2194 G.(CC4 9.001) C.)»C4 


t 1.494t-02 

Vise 1,4J4( 32 

«|J l.S2t( 00 

mi 1.0141 til 

on i.4s*t 01 

II 1.42*1-02 

m -i.iJM-oi 


lUAK 42TI 

1 V I 

10.00 0.00 0.09 

O.CO >9.00 0.00 

0.90 0.00 -9 .09 


,1041 .11494 .0441 

,2)04 .174) .0244 

,214* .C249 .1424 


,|4js .Of)2 .07)2 

.1114 .1944 .9*2* 

.1742 .0424 .1044 


N Rij/mi 

u o.m 

If 0.001 

1) 0.001 


mine. o)-)TRtii evo(»7 


LJ/M 11 
n 9.0424 
12 -.099C 

11 0.000< 
22 0.0)41 

2) C.0004 

)J O.OI4» 


\> 1? 22 24 

.'JCI -.0001 0.1114 0.0604 

,0211 -.0(04 O.OOOl -.0604 

O.0964 -.0211 O.OOll 3.0007 

0.0411 C.0502 0.1474 -.0014 

iC(l2 -.tiSiC 0.0017 -.174) 

O.ftl) 9.9694 0.4411 -.'d:! 


! 1.2011 - :i 

vise 1,4141-02 

mi ).2R4t 

fit 1.2441 01 

Oil 9.7C0I n 

ti *.iiu-:j 

111 -2.041I-C4 


MtAd 0t7C4l«*7IC»' »*7‘ 
Jl T 2 

2Q.C0 0.00 0.00 

0.00 -19.00 0.00 

n.co 0.00 -10.00 


IhTl^aU $(4(11 
* » 2 
.lJU .0*44 .04*4 

.0441 .llte .0414 
,044* .0449 .1)14 


,04)7 .0402 .9104 

,CIC4 .04*8 .044 

.(lit' .04M .09)* 


ij m;/mi 

It 0.244 

12 O.OdO 

13 O.OCt 
22 O.ko) 

21 -g.ooa 

1) 0.M4 


IJ/Mt 11 
11 9.64$) 

U 9.69C2 
) 9.6000 

I 6.2i«; 
.1 9.0000 

n 0.2114 


• I IHOLOS- StM SI BVjOOtY 


- f4t$Sti4t-S14*14 7I4S0* 


.CCCO 0.0904 0.2144 

,0414 - .0004 - ,0002 


0.2144 -.900) 
-,9002 -.9900 
-,Cm01 0.9002 


0.CCC4 .9)7^ -.CmOI 0.9002 
0.06(9 O.OOC9 0.1M7 0.9904 

-.ten -.0«(7 O.CtiOl -.1?)) 
C-.CCtI 0.0C12 0.7)41 -.0007 


7 4. Tin : 

Vise 1.41*1-2 
•II 2.4t*( c 

• M 1.9141 I 

on 4. 1 0*1 \ 

M I S 911-02 

1 1 1 .4,S4)f- " 


• |«N CMC«"47K4 •*’l 
» T I 

IW.Cr J.OO «.9fl 
O.C) -10.00 0.09 

0.0) C.OO -10.90 


tM{w44k SCk4(S 
I 7 2 

,11)4 ..SM4 .09>4 

.22)1 .1*94 .9144 

.2191 .9204 .11*4 


.1)41 .9499 .0*19 

.144* .9497 .3934 

.1441 .9109 .(144 


■ r 4iirLQ$-$7*M9 ItiCOH 


-0.C92 -0.991 

O.On 0.99* 
-O.OOl *0.091 
0.0J4 0.947 


ms • •!$ ' 

.124 O.0$4 

,.C9l 9.091 

’.(Cti -O.wOO 
I.C4* -0.917 


MIS$V«C-17«»lli Tfl«04 


li/Nt 11 
1*. 0.0H2 

11 9.0002 

}1 9.4002 

ii 4.9441 
2) 0.9042 

n 9.94»* 


.C((7 -.6602 1.1819 

-.5110 -.)69l 9.6002 

e.(co) -.0171 e.X(4 

d.0C(7 O.04U 9.2049 

-.{(12 -.9CJ4 O.OOl) 

3. >(10 9.9019 0.4774 


.1689 9.0094 

.0002 -.0092 

.K(4 8.9999 


I *,I99(-0I 

VlK 1.41*1-62 

fill ).09*f 0) 

mi 4.59)1 e< 

on 1.77)7 91 

I I f .t)9r-ei 

in -1.4)41-6) 


MAN DMOONATICfi 8*71 
*82 
10. CO O.CO 0.00 

0.90 ■>)»80 0*90 

0.00 6.00 **1.40 


IN710AAL UAUfS 
A V 2 

.1)92 .0771 .9744 
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APPENDIX B 


DETAILS OF THE NUMERICAL ALGORITHM 


THE EQUATIONS OF MOTION FOR HOMOGENEOUS TURBULENCE 


We begin with the Navier-Stokes equations for a fluid of uniform density and 
viscosity: 


«“i 

— + 


VU 


ijj 


(Bl) 


“i,i “ 0 


The velocity field is decomposed into mean and turbulence components 
ui - Ui(x) + u£(x') , p ■ P(x) + p^(x'”) 

where x' is related linearly to and the mean is restricted to a spatially linear 
form 

Ui * Uijxj (B2) 

where depends only on time* 

The explicit appearance of (but not its gradient) in the equations of 

motion for the turbulence field is eliminated by the use of a spatial coordinate 
system that moves with the mean flow. This eliminates the terms, which contain x 
explicitly, that represent convection of the turbulence by the mean velocity. 

This moving coordinate system is linearly related to an Inertial laboratory 
frame by the relation 


x£ = B.jX. 


(B3) 


where B is determined by the constraint that the new coordinates move with the mean 
flow or~ 


Dx' 3x' 


3x;' 

1 


Dt 


3t 




so that 


k 




ij ' "Ik k,J 

The Navier-Stokes equations then become 


(BA) 
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1 


0 


3u,' 

<Ui,k + + P,1 + -g^ + Uj'Ui^j + Bkj(uiUj)^k + Bkip;k - ^BkJ®£J“i,kt 

Ui,i + “ 0 


Note that differentiation of U is with respect to x while that of is with 
respect to x''. The time variable for both is t. 


To separate the mean from the fluctuation parts we take the ensemble average of 
the Navier-Stokes equations » using the ergodic hypothesis to replace ensemble 
averages of spatially homogeneous terms by their volume averages: 


(Ui.k + Ui.jUj,k)>^k + P,i = 0 

Ui.i = 0 

This gives the equations for the turbulent component* 

3u' 

— + 

BjiUij - 0 


(B5). 


(B6) 


The mean-momentum equation implies that - P must be quadratic in 21 (ignoring 
an arbitrary additive function of time), that is, 


P - Cji^XjXj, 


SO that 


p.i “ - ^ik’^k 


where P is a symmetric tensor* 

The mean-momentum equations then reqv^ire 

’ji.k + ^'l.j^'j.k + Pik = 0 

When the mean deformation is decomposed into strain and vorticity tensors 

= ‘’ij ^"ij 

we find upon separation Into symmetric and antisymmetric parts that 

s + S" + n^ + P“0 

Q + sn + r.s = 0 

This is simply a stateBeol. tjD.it the atrAlP rate S may be specified as an arbitrary 
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function of time subject only to the constraint - 0. The second equation 
determines n, and the first equation determines the mean pressure field. The 
initial value“of fi is arbitrary. 


The basic equations for the turbulence field (B6) can be solved in a 
straightforward manner in principle, but in practice one is usually constrained by 
computer memory size, and it is advantageous, and in our case absolutely necessary, 
to reduce the number of words per node to the minimum. For this reason, the general 
case of (B6) with arbitrary mean-deformation tensor was not implemented. The 
current program allows a velocity deformation tensor of the form (rotation is handled 
as a special case) 


'i.j 



a + b + c 


(B7) 


If the mesh is given at t ■ 0 by 

Bn(0) = , B22(0) = &2 . 

the time variation of the metric tenser 


B33(0) = S3 , rest of B zero 
3 is determined by (B4) as 


- /* *■ adt 
, "0 

cdt 



(B8) 


®12 


-s(0) 

S1S2 


B22(b) 



Bn(t)dt 

, rest of B » zero 


The use of Si other than unity allows the mesh spacing to vary with direction. We 
have made use of (B7) to find that 


s(t) 


fhdt 

s(0)e ^ 


s(0)S3 

b^TuT 


RAPID-DISTORTION THEORY 


Substitution of (B7) Into (B6) gives (taking for the moment s - 0) 
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Ut + au + Bi iPx ■ y 

Vj + bv + B22Py • [ nonlinear + viscous terms 

+ CW + B 33 P 2 ■ / 


(B9) 


Biiu^ + B22Vy + B33WJ5 - 0 

In the rapid-distortion limit the strain rates 
integrable, so that the metric tensor B has a 
straining period, 0 < t < T, and (B9) degenerates 


a, b, and c become Infinite but 
finite discontinuity. During this 
to 


Uf + au + BiiPjj " Bn 


where 


r 


p dc 


with similar equations for v and w . 

This suggests the Introduction of the scaled velocities 

(u',v%w') - (B7|u, B-lv, B-lw) 


(BIO) 


and new dependent variables u + v + 4>y> 

that remain continuous during a rapid distortion. In terms of these new variables 
the rapid-distortion process is simply 


u'(T) = u'(0) - 4 >x(T) 
v'(T) * v'(0) - 4>y(T) 
w'(T) - w'(0) - i'jCT) 


(BID 


where ^ (T) is found by applying the continuity condition (B9) at t - ' 
Bn<fxx + bLv^IO) + B 33 W^( 0 ) 


When the shear rate s is not zero, the 
pressure cannot be completely absorbed into a 
The equations analogous to (B9) are 


metric B varies with time and the 
new dependent variable as it was above. 


Ut + au + B; iPy 
V[ + bv + B22Py + Bi2Px 
Wt + cw + B33r2 
BllUx Bv2Vy + B33W2 


nonlinear + viscous terms 


(B12) 


and we take as new dependent variables 
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X « 


(B13) 


Bi2 


w' + 4>s 


In these varlableii equation (BIO) remains valid for the x and z directions, but 
in the y direction 


Vj + bv + B 22 Py + Bl2Px “ ®22 
where, from (B8) , 


®22'^ <*y 


®12 \ d /®12\ 

B22 ‘’’V ' dt \B22/ 


d_ / ® 

it \B22/ Bi 


_d 

d't 


zim «? 
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Bll 


SC that ^22 ^x 


-f— ) 

^ dt \B22/ 


remains bounded during a rapid-distortion. 


The new dependent variables above are therefore continuous through a rapid 
distortion, and, as before, the state of the velocity field following the distortion 
is found directly from the continuity condition (B12). 


The numerical advantage of the new dependent variables is that the strain rate 
no longer imposes a stability limit on the time— advance step size, and that 
rapid-distortion theory is exactly captured by the simulation. 


THE VISCOUS TERMS 


In the moving coordinate system, the Laplacian becomes 

and the diffusion operator in wave space becomes 

^ ■ F-l(k.t) ^ |F(k,t)u,l 

where F is the integrating factor 

log F - vk,^kj J B^jB^j dl 


We again re-define the dependent variables to include the integrating factor and 
eliminate the diffusion term as <1 numerical stability constraint. The use of the 
integrating factor also permits us to carry fewer words per node for some 
time-advance algorithms (e.g., fourth-order Runge-Kutta) than the explicit treatment 
would allow. 
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In addition to the more or less "classical" transformations above, we have made 
a few simple modifications of Che algorithm that pertain to data management. The 
primary consideration here is the reduction of the number of words per node required 
for the simulation. A secondary consideration Is the cost reduction that results 
from lowered operation counts and lowered I/O operations between fast memory and 
backing store. 

The necessity for thes» modif lections and their actual form is largely dictated 
by the size and architecture of the llliac computer, and the data-base structure 
chosen by the author. These modifications are presented for the basic Navier-Stokes 
operator in Cartesian coordinates; their application to the transformed system of 
equations Is straightforward. 


PSEUDODERIVATIVE 


Wc u«e as dependent variables 


Q(0,k2,k3> , Ikiu(k) , 

vCki.O.k^) , Ik^vCk) , 


w(ki,k2,0) , ik3w(k) , 


ki 0 

k? * 0 
k3 0 


rather than the usual Fourier coefficients: 

u(k) , V(k) , ^(k) 

In physical space the equivalent Information is 


u(0,k2,k3) 



u(x,y ,z)dx 


etc. 


ikju(k) ■» >■ 


3x 


u(x,y,z) 


etc . 



The basic variables are then Uj,, Vy, and w^. The integrals determine the constants 
of integration required to recover u, v, and w by Integration. These variables are 
rather natural for the Navier-Stokes equations because their use leads to 


^xt 

+ 

(uu) 

4- 

iuv)xy 

+ 

( ) X35 

Vyt 
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Uw)xy 

+ 

(VV)yy 

+ 

f 

N 

Wzt 

4- 

(uw) X 2 

4“ 

(V-Wlyj, 

4- 

(ww)zj. 


in which the "matrix" of nonlinear terms is symmetric, even when the differentiation 

A3 









operations are included ^ In the usual variables we would have nine different 
nonlinear terms rather than the six above. 

An additional benefit of these variables is the ease of enforcing the continuity 
condition. Thus it is practical to discard w^ from the database on backing store 
altogether since it can be simply recovered from and Vy. Of course we must 
save S(ki,k2»0), but this is only a two~dimensional array. 


MODIFIED PRESSURE 


The nonlinear terms are formed (in physical space) when (x,z.) planes are in 
fast memory. Thus we can form the nonlinear terms: 

(uu)jjjj + (uw)jj 2 » « (vv) , (vw )2 > (nw)xz + 

The y differentiation must wait for a pass over the database bringing (x-y) planes 
into fast memory, and this prevents further summing of terms. Thus we would have to 
carry five words per node to calculate the nonlinear terms. This can be reduced to 
four words per node by simply introducing the modified pressure p + v . The 
required nonlinear terms are then 

(u^ - + (uw)j^2 , , (vw)^ j , + (^”)xZ 

Because cf Tlliac constraints, an implementation of five words per node for the 
calculation of the convective terms would have been difficult and considerably less 
efficient than that for four words per node. 

The disadvantage of the modified pressure is that we have lost the ability to 
form output statistics containing pressure, because we only have the modified 
pressure and not the • v‘ information necessary to recover physical pressure from it. 
When the pressure is required for output statistics a separate calculation must be 
made . 


FINITE FOURIER TRANSFORMS 


The use of fast-Fourier-transf orm techniques in the simulation of isotropic 
turbulence was pioneered by Orszag and Patterson (1972). Because we have departed a 
little from their treatment, it is. worthwhile to cover some of the details of our 
implementation. 

For purposes of analysis it is advantageous to work with general complex 
tiansforms, even though the physical data are real. In practice, of course, the 
complex transform is not used as such because it doubles the time and space 
requirements. (One can, however, use the complex algorithm on two real fields 
simultaneously, or use a half-length complex transform on full-length real data.) 


The only subtle aspect of the use of discrete Fourier transforms Is the aliasing 
problem. Although it has been discussed by many other authors, it seems appropriate 
to give a short exposition of it here. 


The one-dimensional discrete Fourier transform of length M is defined as 


Ma 


n 


^ -2iTi(jn/H) 

L aj e 

j-0 


with inverse 


M-l 

E 


n»0 




^2TTi(jn/M) 


Suppose we have data aj given by 


a 


j 


^ 2irl(j/M)k 
p o 


a wave having integer wave number k.. The discrete transform of the.ie data gives for 
M values of n ; 


= M6<St(k - n)mod M] 
j“0 


SO that 


= 6iS[(k - n)mod Ml 


Thus a wave having wave number k in physical space Is transformed into the 
Fourier coefficient at wave number k mod M in wave space. This can also be seen by 
noting that the discrete data 




B e 


2rri(j/M)k 


e 


:7ii(j/M)(k+NM) 


for any integer N. 

Thus the data at 
same for any value 
number k • 


M equidistant points 
of N. Wave numbers 


of the period In 
k+NM are said to 


physical space Is the 
be "aliased" to wave 


The aliasing problem occurs when operations on data (usually 
physical space) with a wave--number span of M produce a wave-number 

than M. 


occurring in 
span greater 


In the Navler-Stokes equations, the nonlinear (bi-llnear products) terms act to 
Increase the span of wave-numbers. Consider a product in one dimension. 
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ijbj 


M-1 

£ 

n“0 

M-1 

£ 

n=0 


2iri(jn/M) c 2iTi(jm/M) 

3jj e 2-r ® 

m=0 

. c 2 iri(j/M)(n+m) 

itpO 


which has the discrete Fourier transform 

M- 1 M- 1 , 

2k - £ £ anb„6((n+m-k)modM] 

n=0 m~0 


<“ A/r - - ,=r.o":r, 

2-M < k < M. This is illustrated by the ^ f oduce the aliased wave 
product of waves with positive wave-numbers similar way products of 

at n-hn-M rather than at the correct wave t^a positive wave-number in 

waves both having negative wave-numbers can be aliased to a posi 

the result. 



The 

different 


alias error resulting from bi-linear 
methods (Patterson and Orszag (1971)). 


products 


can be 


removed by two 


ALIAS REMOVAL BY TRUNCATION 


we want to form the length M /esuU°"In'a Fouiter^serle^o? 
M (l-M/2 < k < M/2). As ^ allL-free product can be obtained by using 
length 2M-1. It is then obv.ous ^hat ircarried out as follows: (D form 
a Fourier transform of length _2M. l_M/2 < k < M/2 and 

zeroes for the remaining 


kl (2) irv^r? tie transform to obtain the length 2M 










physical space; (3) form the length 2M product ab ; end (4) return to length 2M 
wave space » using a length 2M transform* At this point we have the length 2M 
alias-free Fourier transform of the product. We now simply discard the modes in the 
result outside our length M original data .to. obtain the length M aliaa-free 
product . 

The same result can be achieved using transforms of length 3M/2 rather than 
2M, with a resulting time and space savings* Because we are going to retain in the 
results only those waves l-M/2 <k<H/2 we pose the following problem; What is the 
shortest length (L) transform that will not alias the product into modes 
l-M/2 £ k. £ M/2 ? The limiting product is that between modes both having wave-number 
M/2 resulting in wave-number M aliased by the length L transform to wave-number M-L> 
which must be outside the span l-M/2 ^ k ^ M/2* Thus we have 

M-L < l-M/2 or L > 3M/2-1 

Similarly I the lowest wave-number I-M/2 at worst produces wave-number 2-M in the 
result and is thus aliased to L+2-M which must be greater than M/2. This requires 
L > 3M/2-2. 

In practice* the length of the Fourier transform (M) is determined by computer 
size and the wave-numbers carried are reduced from l-M/2 ^ h £ M/2 to -M/3 k £ M/3 
in order to achieve alias- free products* C^ote that by M/3 we mean the truncated 
integer value.) This procedure is popularly referred to as the 2/3 rule. 


ALIAS REMOVAL BY PHASE SHIFTS 


We can rewrite the transform of an aliased product as 


n,m n*m 

n+m=k n+nv=k±M 


where the first sum is the alias-free result* and the second sum is the alias error 
at wave-number k. Consider the effect of a shifted physical space grid. In wave 
space the shift is achieved by multiplying each Fourier mode a^^ by the phase factor 
ginA^ If we form the product on the shifted mesh and then shift the result back 
to the original mesh we obtain 


or 



inA.^ imA . 
e b„, e -h 


E 

n+m=k iN 


inAc 
e bm 



Ck 


n+tn=k n+m=kiM 
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The alias-free sum is unchanged by the shifts, but the alias f 
hv the Dhase factor This offers several possibilities for eliminating the 

-.r: :;n xj-- 


±iMA^ 

= -e 


and the alias-free sum is one half the sum of the two evaluations. 

ThU Mthod »Uovs o„e to retain all of the Fourier mode, aeeoolated »lth . 
length n transform, but requires two evaluations. 

ret.ulre'rouUe’’"™" “to;, .lelhe'^e f “ 

Random value of taken from a uniform distribution around the unit circle. 

The current code, as we will see shortly, uses both truncation and shifts in 
concert to effectively eliminate alias error at little cost. 


MULTIDIMENSIONAl. ALIAS ERROR 


The -,Uiai»ene.o„el discrete Fourier thrillesed Se^e-a^rr'ire 

:i:::ir;:ctiL th ueve-uu.hete 

L’’jrerL.l:?n/re“co;erLd.^Thetle,the result he 

aliased nr not In each of. the three directions independently. 

When independent shifts are applied to each direction, we have 


Sc .0 


OiSioO 


+ A ,S, 


) 1 xO 


+ 1 ‘ 




00 1 




+ 02t\,Soii + ei0 3 Sioi + 0162038111 


where 


A 8 


... I 




A. 


®1 ■ (etc. for 62 » 03 ) are the phase factors, 


Sooo ■ 2 a(n) 6 (m) 

m+n*k 


'100 “ a(n)b(m) , 

i!^=k+Miei 


is a unit vector in the x direction 
etc. for Soio, Sqoi 


'110 “ 3(n)b(m) , g 2 is a unit vector in the y direction 

i_i+n-k+Miei+M 2 e 2 etc. for Eji,, Siqi 


®111 ■ a(n)S(tn) 

m+n^*k+M 


■> « alias-free result is Sqoo; Siqo. Sqio* and Sqoi are sums aliased in the 1 , 

7« f ** 3 directions, as indicated by the subscripts; S^o. Squ, and Sjoi are aliased 

in two directions, and Siij is aliased In all three directions. 

We use a combination of truncation and phase shifting similar to that of 
Orszag CO eliminate all of the alias errors. Recall that truncation by 
the ^3 rule consists of zeroing wave-numbers ki with |ki|> M/3 prior to inversion of 
the transform, and their discard upon return to wave space following formation of the 
bi -linear product . 

dimensions truncation of wave vectors having any component greater 

han M/3 results in the elimination of all the alias-error sums S for the remaining 

wave ve._Lors. If only those wave vectors having two or more components greater than 
M/3 are truncated, only the double and triple aliases are eliminated. If only those 
vectors having all three components greater than M/3 are truncated, only the triple 
alias sums are eliminated. We have followed Patterson and Orszag and eliminated the 
double and triple aliases by truncating wave vectors having two or more components 
greater than M/3 (Patterson and Orszag truncate modes with > 2fM/3)^ a more 
severe truncation!. We are then left with only the single aliases: 

^k " ^000 ®lS;oo ■h 62S010 *^38001 

These can be eliminated by summing evaluations from shifted grids as outlined 
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previously. If the exact elimination of the single-alias error is not required, we 
can avoid the extra evaluation and perform the phase shifting as the time-advance 
proceeds. Thus the alias error from one step will be nearly cancelled at the next 
step. In the Runge-Kutta methods used in the current algorithm, the residual alias 
error is of the order of the square of the time-step (h). Each cancelling pair of 
evaluatious is made at the same time level of the advance, and the shifts in the 
three directions are uncorrelated with one another, that is 

e(i) = -0^^) 

where are three random complex numbers chosen at each time-level from a 
uniform distribution around the unit circle. The alias errors act as a random 
forcing function of small amplitude (which is formally O(h^), but also dependfs on the 
energy spectrum). 


RE-MESHING THE FIELD 


The computational coordinates move with the mean flow and can become greatly 
distorted Irom their initial form. When the mean flow is an irrotational strain, the 
coordinates remair orthogonal, but the mesh aspect ratio grows Indefinitely. Because 
they are influenced by the nonlinear terms in the Navier-Stokes equations, the length 
scales of the turbulence field are not simply stretched and rotated by the mean 
field. As the mesh becomes distorted, the range of scales of the turbulence in some 
direction may not be contained within the resolved range of scales in that direction 
even though ir. was at earlier times. We have not yet studied the general re-mesh 
problem, and the current code re-meshes only the shear flows. The process is 
illustrated by the following diagr^'m for the case gj The extension to 
general 3 is straightforward. 



The calculation begins at st • 0 on a Cartesian mesh, shown above as the solid 
vertical line. At time st 1/2 the mesh has become skewed to the right as shown by 
the dashed line. The field is then interpolated onto the mesh skewed to the left 
(the dotted line), and the time-advance proceeds. At time at * 1, the mesh is again 
orthogonal, and at st'3/2 the remesh process is carried out again. 
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K '' r \ / V w. * A 


The re"-meBh process in physical space appears simple , but In fact it results In 
alias errors unless those errors are explicitly eliminated. Denoting the coordinates 
before and after re-mesh by superscripts (1) and (2), we define the re-mesh process 
as 

u(2)[x^C2)(x)] « u(0[x^(0(x)] 


WUCJLC , . > V 

+ x£ ej 

£ is the X unit vector 

so that 


- X2 fill 


To carry, out the reuiesh in physical space we would sltnply shift the data in X“Z 
planes right (end-around due to the periodicity of u) by the amount y. The first 
plane is not shifted; the second plane is shifted by one node, etc* The (M+1) 
plane, which is implicitly equal to the first plane by periodicity in y and not 
actually present, would be shifted by M nodes (one full period) and remain equal to 
the first y plane* In terms of Fourier modes the re~mesh operation is 




Thus a wave with wave-number (ki,k 2 ,k 3 ) on the old mesh has wave-number 
(kj ,k 2 +k; ,k 3 ) on the new mesh and is aliased if k 2 "i‘ki > M/2* 

The alias errors during re-mesh are eliminated by inverting the length M wave 
space data u to a length 2M physical space in the y direction, leaving the x and z 
directions in an M length wave space* The data are then shifted in x, using the 
appropriate phase factor, and the shifted result is returned to wave space in y using 
2M length transforms* The data for wave-numbers k 2 > M/2 are then discarded and the 
re-mesh is complete* Because we have discarded modes, we lose both energy and 
dissipation (both being positive-definite), the loss of dissipation being by far the 
larger because we lose primarily high— "rfave— number information. 


THE INITIAL CONDITIONS 


The Fourier transform of the velocity field Is initialized to an Isotropic 
state, satisfying continuity, and having a given energy spectrum* 

The continuity condition is satisfied by taking 


SI 




where is the computation vector basis and e£ is any vector basis having 63 
parallel to k. The complex components a and 6 in general are random in amplitude and 
phase, subject only to the constraint that 

J<u«u’^> dA(k) « E(k) =J<aa* + dA(k) 

where E(k) is the desired energy spectrum and the brackets < > denote the expected 
value. We have not implemented this general form because we feared that the 
deviation of the energy spectrum from its expected value could be large at low 
wave-numbers where the sample of modes in A(k) is small. 


The form implemented is 


a 



cos <p , 


^ " V^7Tk2; 


where 0^, 62, and ^ are uniformly distributed random numbers on the Interval ( 0 , 27 r). 
This form results from the stronger constraint on (a, 3) 


E(k) = (aa* + B3*) JdA(k) 

where the energy of each mode separately is required to have exactly the expected 
value. The modes are still random with respect to spatial phase ( 01 , 62 ) and velocity 
component distribution ( 45 )- 

To complete the formulation we need only to relate the basis e£ to the 
computational basis Any basis subject to the constraint 

kc3 * kiei + k2£2 + k3g3 = k 

will do, since rotations about 63 can be absorbed into the random phase We 

choose arbitrarily a basi'- having 

• ej = 0 

which leads to a solution for 
(k[" + 

k(k;^ -f k2*“) ^ £2 ” + k 2 k 362 “ + ^2^)63 


Thus we find 
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(akk2 4 6kik3 \ / |3k2k3 - akkj \ /sCk^ 


2 4 k2^)^^^ 


)- 


The initialization process above does not produce the higher velocity moments 
(e.g., skewness) characteristic of real turbulence. To relax the flow to a more 
realistic state we usually allowed it to decay isotropically until the energy cascade 
was developed and the spectra became fairly smooth. The resulting isotropic fields 
were stored to be used as the starting conditions for various anisotropic runs. This 
procedure was used for all of the irrotational-strain runs in which the experiments 
presumably also had fairly well-developed and nearly isotropic turbulence prior to 
passage through the straining section of the tunnel. For the shear cases, in which 
we were primarily interested in the later stages of development having growing 
energy, we simply ran from a square-pulse spectrum. 
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Figvire 4.- History of the component cnergie.'^ in homogeneous turbulence 
subjected to axisymmetric strain. 
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Figure 17.- Dependence of shear-stress correUtlon on Reynolds number 
in homogeneous shear turbulc.-^ce. 
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Figure 21.- Development of velocity anisotropy of homogeneous turbulence 

in uniform rotation. 
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Figure 22.- Velocity gradient anisotropy of homogeneous turbulence 

in uniform rotation. 
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Figure 28 .“ Comparison of measured and modeled terms in the Reynolds stress 
equations for homogeneous shear turbulence. 
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Figure 29.- Variation of Rotta coefficient with Reynolds number 
in homogeneous shea^* turbulence. 
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Figure 30,- Variation of Rotta coefficient with anisotropy 
in hom^geneoas shear turbulence. 
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Figure 32.- Compa.iaun of measured .and modeled values of the "f ast-pressure 
strain" tensor using the l.avnder-Reece-Rodl model in homogeneoua 

shear turbulence# 90 
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Figure 33.- Variation of the Launder-Reice-Rodi coefficient with Reynolds 
number in homogeneous shear turbulence. 
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Figure 3^^,- Variation of the Launder-Reece-Rodi coefficient with anisotropy 
level in homogeneous shear turbulence, 
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